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INTRODUCTION 

In  many  countermeasure  applications,  it  is  necessary  to  broadcast  broadband  white  noise  to 
influence  a  device  that  has  a  relatively  narrowband  processor  of  unknown  center  frequency.  Due 
to  practical  limitations,  the  excursions  of  the  power  supply  generating  the  broadband  noise  are  hard 
liimted  by  a  peak  factor.  An  important  operational  question  then  is:  to  what  extent  does  this  limit 
the  excursions  of  a  narrowband-filtered  output  of  the  power  supply? 

In  this  note,  the  broadband  white  noise  is  modelled  as  a  finite  sum  of  discrete-frequency 
sinusoids,  with  Gaussian-distributed  amplitudes  derived  from  a  uniform  distribution.  The  resulting 
time-domain  signal  then  is  hard  limited,  and  the  discrete  Fourier  transform  is  applied  (exact  for  this 
type  of  signal)  and  executed  via  the  fast  Fourier  transform.  The  resulting  transform  amplitudes 
then  are  hard  filtered  for  a  particular  choice  of  bandwidth,  and  the  inverse  discrete  Fourier  transform 
is  applied.  Time  samples  of  the  resulting  narrowband  time  signal  then  are  used  to  update  a  histogram 
of  v^ues,  the  amplitudes  of  the  broadband  time  signal  then  are  re-seeded,  and  the  process  repeated 
until  a  sufficient  number  of  signal  samples  are  collected. 

The  resulting  histogram  is  used  to  constmct  a  least-squares  fit  to  a  Gaussian  distribution,  and 
the  resulting  variance  is  compared  to  that  for  the  unclipped  case. 

THE  SIGNAL  MODEL 

The  broadband  signal  has  the  form 

NI2 

B  (n  At)  =  X  [a„  cos(2jinm  At  A/) + b„  sin(27i«m  A/A/)]  (1) 

m  =  1 

where  A/  =  1/7  and  A/  =  1/2F  with  7  the  period  of  the  lowest  frequency  sinusoid,  and  F  the  highest 
frequency  present.  We  will  choose  /V  to  be  a  power  of  2,  and  note  that  I  ^  N/2,  and 

A/A/=  VN. 

The  coefficients  are  chosen  using  the  following  algorithm: 

U, eU[0,l]  ,  Ujemi]  (2) 

V,  =  2I/,-l  ,  V^  =  lU^-\  (3) 

that  is,  Ux  and  I/j  are  chosen  from  a  distribution  uniform  on  the  interval  [0, 1],  and  are  used  to 
generate  a  distribution  uniform  on  [-1, 1],  of  which  V,  and  are  samples(most  computers  generate 
the  distribution  uniform  on  [0, 1]  by  means  of  a  pseudo-random  number  generator). 

Define  S  =  V?  Vj .  If  S  >  1  then  select  new  samples  C/j  and  Ui  Trom  the  uniform  distribution 
and  repeat  the  process.  The  resulting  variable  S  is  distributed  uniiormly  on  [0, 1].  Next,  define 
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and 


*.=v,V^ 


(4) 


(5) 


The  resulting  variables  andX2  have  Gaussian  distributions  with  zero  mean  and  variance  a. 

For  a  given  broadband  signal  sample,  a  pseudo-random  sequence  is  seeded,  and  the 
coefficients  a„  and  b„  are  selected  using  the  above  procedure  for  1  <  m  ^  N/2.  When  a  new  signal 
sample  is  constructed,  the  sequence  is  re-seeded. 

When  a  peak  factor  P  is  specified,  the  clipped  signal  has  the  form 

CinAt)=B(n6i)  <aP  ,  C(nAr)  =  aP  otherwise.  (6) 

The  discrete  Fourier  transform  of  the  clipped  signal  is  given  by 

D(/A/)=  i  (7) 


with  the  inverse  transform  given  by 


C(nA^)  =  ^^D(/A/)e^'"'^ 
N  t=i 


(8) 


The  validity  of  this  transform  pair  can  be  established  by  means  of  the  identity 

N 


I  e 


2Ki{n-k)mfN  _ 


=  S  5, 


'n,k+IN- 


(9) 


NARROWBAND  STATISTICS 

The  discrete  Fourier  transform  and  inverse  transform  can  be  executed  via  the  fast  Fourier 
transform  algorithm.  A  subroutine  for  performing  this  algorithm  is  included  as  part  of  the  code 
given  in  Appendix  B.  The  code  is  written  in  Hewlett-Packard’s  Rocky  Mountain  Basic  Version 
3.0.  The  main  body  of  code  constructs  a  histogram  for  a  particular  bandwidth  and  peak  factor 
choice.  For  a  particular  broadband  signal  sample,  the  time  signal  is  transformed,  and  then  hard 
filtered  by  elii^ating  the  transform  coefficients  outside  the  chosen  band.  In  ^ing  tlus,  it  is 
important  to  observe  the  redundancy  and  symmetry  of  the  transform  coefficients.  In  particular, 
because  we  are  using  a  real  time  series  as  input,  the  transform  coefficients  satisfy 

Re{DilAf)}=Re{D(lN-l]Af)}  ,  /  =  1 . N/2  (10) 
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and 


/m{D(/A/)}  =  -/m{D([A/-/]A/)}  ,  l  =  l,...,NI2.  (11) 

Consequently,  a  set  of  bandpass  indices  is  selected  symmetrically  about  /  =  N/2.  Once  the  filtering 
is  accomplished,  the  coefficients  are  inverse  transformed  to  produce  the  associated  time  signal. 
Each  time  sample  is  examined  and  assigned  to  the  appropriate  histogram  bin.  Once  this  is  done, 
the  process  is  repeated  until  a  sufficient  number  of  samples  is  accumulated. 

For  explicit  analysis,  we  have  chosen  F  =  2048,  and  A/= 4,  so  N  =  1024,  T  =  1/4  and 
At  =  1/4096.  The  broadband  signal  is  chosen  to  have  o  =  1,  so  the  amplitudes  a„  and  b„  are  chosen 

via  the  above  procedure  using  a  c  of  Vl/512.  For  the  filter  we  consider  bandwidth  ratios  of  1:40, 
1:4, 1:2, 3:4, 7:8, 125:128,  and  1:1,  centered  on  the  frequency  1024.  The  histogram  is  constructed 
as  an  array  of  2M  elements  with  element  M  corresponding  to  zero  amplitude.  In  some  cases,  we 
have  used  M  =  100  and  in  some  cases  M  =  1000.  It  doesn’t  appear  to  make  much  difference  in  the 
resulting  least-squares  Gaussian  fit  The  matrix  index  change  corresponding  to  a  for  the  unclipped 
case  is  arbitrarily  chosen  to  be  Af/ 10  and  this  number  is  multiplied  by  the  signal  value  and  divided 
by  the  square  root  of  the  bandwidth  ratio  and  the  result  added  to  M  to  give  the  appropriate  histogram 
index.  This  will  give  the  same  distribution  for  the  undipped  case  for  all  bandwidth  ratios.  Then 
the  effects  of  the  peak  factor  will  appear  as  departures  fit>m  a  Gaussian  distribution,  and 
modifications  of  the  heights  and  widths  of  the  Gaussian  fits. 

The  distributions  and  Gaussian  fits  for  a  bandwidth  ratio  of  1 :40,  for  peak  factors  of  1 .7, 1 .5, 
1.25, 1.0,  and  0.5  are  shown  in  Figure  1.  Here,  and  in  Figure  2,  the  abscissa  is  labelled  in  units  of 
o  for  the  undipped  case.  The  smooth  curves  are  produced  by  the  least-squares  fitting  procedure 
described  in  Appendix  A.  It  is  clear  from  the  comparisons  given  in  1(b)- 1(0  that  the  distributions 
are  Gaussian.  The  superposition  of  fits  in  1(a)  shows  that  tiie  effect  of  increasing  the  peak  factor 
is  to  reduce  the  variance  and  increase  the  distribution’s  peak  amplitude.  It  is  clear  from  the  curves 
in  1  (a)  that  imposing  a  peak  factor  of  1 .7 ,  or  even  1 .5  has  very  little  practical  effect  on  the  resulting 
narrowband  distribution.  Note  that  for  peak  factors  of  1.5, 1.25,  and  1.0,  we  have  used  M  =  1000, 
and  have  203  samples,  whereas  for  peak  factors  of  1.7  and  0.5,  we  use  M  =  100,  and  have  146 
samples. 

To  determine  when  the  peak  factor  does  have  a  significant  impact  on  the  shape  of  the 
distribution,  we  start  at  the  other  extreme,  the  distribution  for  the  full  bandwidth  clipped  signal, 
and  reduce  the  bandwidth  in  slight  increments.  In  Figure  2,  we  show  the  distributions,  for  peak 
factor  1.7,  for  bandwidth  ratios  1:1, 125:128, 7:8, 3:4, 1:2,  and  1:4.  In  the  1:1  case  shown  in  2(a), 
as  expected,  the  distribution  has  a  delta  function  behavior  at  ±1.7a.  In  the  next  case,  that  for  the 
slight  reduction  to  a  125: 128  ratio,  shown  in  2(b),  the  distribution  is  still  strongly  multimodal.  The 
distribution  is  still  noticeably  non-Gaussian  f^or  subsequent  bandwidth  r^uctions,  shown  in 
2(c)-2(e),  but  for  a  ratio  of  1:4,  shown  in  2(f),  the  distribution  is  nearly  Gaussian  with  a  variance 
of  92  percent  of  that  of  the  undipped  case.  In  these  examples,  the  number  of  samples  varies  as 
follows:  1:1(303),  125:128(177),  7:8(178),  3:4(178),  1:2(223),  and  1:4(123). 
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ed  after  the  Peak  Factors  are  Given  Relative  to  an  Un- 
ped  Value  of  1.  (a)  Gaussian  Fits  Compared,  Including 
Unclipped  Case.  (b)  1.7,  96.8^.  (c)  1.5,  93.1^.  (d) 
.  84.9*/..  (e)  1.0,  75.8^.  (f)  0.5,  44.5*/i. 
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5 


1.7,  for  Various  Bandwidth  Ratios:  (a)  1:1,  Cb)  125:120, 

(c)  7:8,  (d)  3:4,  (e)  1:2,  (f)  1:4.  The  Distribution  is 

Strongly  Non-Gausslan  Initially,  but  Approaches  a  Gauss  1  an 
Distribution  of  Reduced  Variance  as  Bandwidth  Decreases. 


APPENDIX  A 


LEAST-SQUARES  FITTING  OF  A  GAUSSIAN 
DISTRIBUTION  TO  A  HISTOGRAM 
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Given  a  series  of  variable  and  probability  values  {j:;,  y,}  where  /  =  1, . . A^,  we  would  like 
10  establish  the  parameter  a  for  the  representation 


1 


^|2K<5 


(Al) 


By  taking  the  natural  logarithm  of  both  sides  of  (Eq.  Al ),  we  can  produce  the  equivalent 
equation 


M;  =V.W-t-lnH' 


(A2) 


where 


M,  =  21n(V^y,) 


V.  =  -X: 


and 


(A3) 

iA4) 


(A!) 


The  least-square  fit  is  accomplished  by  numerically  solving 


(w,  -  v.w  -  In  w)*  =  0 
aw 


{A6) 


or,  expanding,  and  incorporating  all  of  the  data  by  averaging. 


2  ,  <Ui>  Inw 

<  V.  >  W  -K  V;  >  Invv - -I- -  =<  M:V;  >  -  <  V:  >  . 

w  w 


(AS) 


A  subroutine  for  performing  this  process  is  given  on  the  following  page. 
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Rocky  Mountain  Bahic  3.0  Code 


10  SUB  Gaushisfit(X(*),Y(*),Npts31>S2,Nmax3) 

20  ! 

30  IGaussian  fit  to  y(x)  using  Npts  points  to  produce  variance  S 
40  !S1  and  S2  bound  the  search  interval  for  S,  Nmax  is  the  number 

50  !of  partition  refinements  allowed. 

60  ! 

70  OPTION  BASE  1 
80  DIM  U(2000).Vv2000) 

90  REDIM  U(Npts),V(Npts) 

100  Lc=LOG(2*PD 

110  Nav=Npts 

120  FOR  1=1  TO  Npts 

130  IFABS(Y(I))<1.E-10THEN 

140  ua)=o 

150  va)=o 

160  Nav=Nav-l 

170  ELSE 

180  U(I)=Lc+2*LOG(Y(I)) 

190  Va)=-X(I)*X(I) 

200  END  IF 
210  NEXT  I 
220  A=DOT(V,V)/Nav 
230  B=SUM(V)/Nav 
240  C=SUM(U)/Nav 
250  D=DOT(U,V)/Nav-B 
260  Xl=l/S2/S2 
270  Xu=l/Sl/Sl 
280  N=0 

290  YO=A*Xl+B*LOG(Xl)-C/Xl-(-LOG(Xl)/Xl-D 

300  Dx=(Xu-Xl)/10 

310  FOR  Xx=Xl  TO  Xu+Dx/2  STEP  Dx 

320  Yy=A*Xx-t-B*LOG(Xx)-C/Xx+LOG(Xx)/Xx-D 

330  IF  Yy*Y0<0  THEN 

340  N=N-i-l 

350  IF  N<Nmax  THEN 

360  Xl=Xx-Dx 

370  Xu=Xx 

380  GOTO  290 

390  ELSE 

400  S=l/SQR(Xx) 

410  GOTO  450 

420  END  IF 

430  END  IF 
440  NEXT  Xx 
450  SUBEND 


A-2 
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APPENDIX  B 

ROCKY  MOUNTAIN  BASIC  3.0  CODE  FOR  SIGNAL  SAMPLE 
GENERATION,  FILTERING,  AND  HISTOGRAM  FORMATION 
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10  OPTION  BASE  1 

21  DIM  Xl(512),X2(512),C(512),S(512).Br(1024),Bi(1024)JI(2000) 

30  MATH=(0) 

40  Samp=0 

50  INPUT  "INPUT  PEAK  FACTOR", P 

5 1  !Peak  value  is  P  times  the  broadband  rms  value 

60  INPUT  "FREQUENCY  INDEX  FOR  NARROW  BAND  FILTER",Ifil 

61  IMidpoint  of  the  digital  filter,  In  units  of  4  hertz 

70  INPUT  "READ  IN  A  PRIOR  fflSTOGRAM(Y/N)?",Ph$ 

80  IFPh$="Y"THEN 

90  INPUT  "MSUS  OF  INPUT  FILE?",Mh$ 

100  INPUT  "FILE  NAME?",H$ 

1 10  ASSIGN  @Filein  TO  H$&Mh$ 

120  ENTER  @Filein;P,Samp,H(*) 

130  ASSIGN  (SFilein  TO  * 

140  INPUT  "DOES  FILE  EXIST  ON  MSI(Y/N)?",Fe$ 

150  IF  Fe$="Y"  THEN 
160  GOTO  240 
170  ELSE 
180  GOTO  230 
190  END  IF 
200  ELSE 

210  INPUT  "NAME  OF  OUTPUT  FILE  FOR  HISTOGRAM?  ",H$ 

220  END  IF 

230  CREATE  BDAT  H$,2002,8 

240  ASSIGN  @File  TO  H$ 

241  INPUT  "INPUT  THE  HALF  WIDTH  OF  THE  FILTER  IN  UNITS  OF  4  hertz"  ,Dlf 

242  INumber  of  lines  retained  on  each  side  of  the  midpoint,  Ifil 

243  Sig2=l/512 

250  ISigma  squared  for  4  hertz  interval  if  1  for  2048  hertz 

251  RANDOMIZE 

260  !  Re-seed  for  each  broadband  signal  sample 

270  FOR  1=1  TO  512 

280  U1=RND 

290  U2=RND 

300  V1=2*U1-1 

310  V2=2*U2-1 

320  S0=V1*V1+V2*V2 

330  IF  S0>=  1  THEN  GOTO  280 

340  Rt=SQR(-2*Sig2*LOG(SO)/SO) 

350  Xia)=Vl*Rt 
360  X2(I)=V2*Rt 

370  NEXT  I 
380  Samp=Samp+l 
390  PRINT  Samp 
400  Timesig:  ! 

410  FOR  It=l  TO  1024 

411  T=(It-l)/4096 

420  INyquist  spacing  for  2048  maximum  frequency 

430  FOR  1=1  TO  512 

440  C(I)=COS(2*PI*4*I*T) 

450  S(I)=SIN(2*PI*4*I*T) 

45 1  !Line  spacing  is  4  hertz 
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460  NEXT  I 
470  Bi(lt)=0 

480  Br(It)=DOT(Xl,C)+DOT(X2,S) 

490  IF  ABS(Br(It))>P  THEN 

500  Br(It)=SGN(Br(It))*P 

501  !  Cli p  the  signal  sample 
510  END  IF 

520  NEXT  It 

530  CALL  Fft(Br(*),Bi(*),1024, 10,-1) 

531  ! Get  discrete  transform  coefficients 
540  MAT  Br=  (l/1024)*Br 

550  MATBi=(l/1024)*Bi 
560  FOR  1=1  TO  513-Ifil-Dlf-l 
570  Br(I)=0 

580  Bi(I)=0 

590  NEXT  I 

600  FOR  I=513-Ifil+Dlf+1  TO  513+mi-Dlf-l 
610  Br(I)=0 

620  Bia)=0 
630  NEXT  I 

640  FOR  I=513+Ifil+Dlf+1  TO  1024 
650  Br(I)=0 

660  Bi(I)=0 

670  NEXT  I 

671  IBandpass  Filter 

680  CALL  Fft(Bre),Bi(*), 1024, 10,-1) 

68 1  ICreate  the  set  of  time  values  for  this  narrowband  signal  sample 
690  FORl=lTO1024 

700  Ix=INT(Br(I)*100*SQR(2048/(2*4*Dlf))+l000) 

701  INormalize  such  that  the  histogram  is  independent  of  the  bandwidth 

702  !for  the  undipped  case 
710  IFIx>0THEN 

720  IFIx<2001THEN 

730  H(Ix)=Hax)+l 

731  lUpgradethe  histogram 
740  END  IF 

750  END  IF 
760  NEXT  I 

770  OUTPUT  @File;P,Samp,H(*) 

771  lUpgrade  the  stored  histogram 

780  GOTO  240 

781  !Go  back  and  re-seed  and  do  it  all  again 
790  END 

800  SUB  Fft(Xreal(*),Ximag(*)JNTEGER  NJMuOsign) 

810  INTEGER  N2,Nul,K,Kl  W,Kln24* 

820  N2=N  DIV  2 
830  Nul=Nu-l 
840  K=0 

850  FOR  L=1  TO  Nu 
860  FOR  1=1  TO  N2 

870  CALL  Bitrcv(P,K  DIV  INT((2''Nul)+.00I),Nu) 

880  Arg=2*PI*P/N 

890  C=COS(Arg) 
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900  S=SIN(Arg) 

910  K1=K+1 

920  Kln2=Kl+N2 

930  Treal=Xreal(Kln2)*C-HXimag(Kln2)*S 

940  Timag=Ximag(Kln2)*C-Xrcal(Kln2)*S 

950  Xreal(Kln2)=Xreal(Kl)-Trcal 

960  Ximag(Kln2)=Ximag(Kl)-Timag 

970  Xreal(Kl)=Xreal(Kl)+Treal 

980  Ximag(Kl)=Ximag(Kl)+Timag 

990  K=K+1 

1000  NEXT  I 

1010  K=K+N2 

1020  IF  K<N  THEN  GOTO  860 

1030  K=0 

1040  Nul=Nul-l 

1050  N2=N2  DIV  2 

1060  NEXTL 

1070  FORK=lTON 

1080  CALL  Bitrev(I,K- 1  ,Nu) 

1090  1=1+1 

1100  IFI<=KTHEN 

1110  Trcal=Xreal(K) 

1120  Timag=Ximag(K) 

1130  Xreal(K)=Xrcal(I) 

1140  Ximag(K)=Ximaga) 

1150  Xreal(D=Trcal 

1160  Ximag(I)=Timag 

1170  END  IF 

1180  NEXTK 

1190  IF  Isign=l  THEN 

1200  FOR  K=1  TO  N 

1210  Ximag(K)=-Ximag(K) 

1220  NEXT  K 
1230  END  IF 
1240  SUBEND 

1250  SUB  Bitrev(INTEGER  IbitrJ,Nu) 

1260  INTEGER  U1,J2 

1270  J1=J 

1280  lbitr=0 

1290  FOR  1=1  TO  Nu 

1300  J2=J1  DIV  2 

1310  Ibitr=Ibitr*2+(Jl-2*J2) 

1320  J1=J2 

1330  NEXT  I 
1340  SUBEND 
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